Sampling from the Bayesian Generator

We now discuss how to use the BayesianGenerator to propogate uncertainty.

As usual we start by bringing in packages

using Random, MarkovChainHammer, Statistics, LinearAlgebra
using MarkovChainHammer.BayesianMatrix
using MarkovChainHammer.TransitionMatrix: generator
using MarkovChainHammer.Trajectory: generate
Random.seed!(1234)

Random.TaskLocalRNG()

and starting off with a generator that generates a stochastic process

Q = [-1.0 4/3 2; 1/4 -2.0 1; 3/4 2/3 -3.0]
dt = 0.01
markov_chain = generate(Q, 10000; dt = dt)'

1×10000 adjoint(::Vector{Int64}) with eltype Int64:
 1  1  1  1  1  1  1  1  1  1  1  1  1  …  1  1  1  1  1  1  1  1  1  1  1  1

We construct the BayesianGenerator object from a prior distribution

number_of_states = 3
prior = GeneratorParameterDistributions(number_of_states)
Q_bayes = BayesianGenerator(markov_chain, prior; dt = dt)

BayesianGenerator 
prior variance  
3×3 Matrix{Float64}:
 1.0       0.416667  0.416667
 0.416667  1.0       0.416667
 0.416667  0.416667  1.0prior mean  
3×3 Matrix{Float64}:
 -1.0   0.5   0.5
  0.5  -1.0   0.5
  0.5   0.5  -1.0posterior variance  
3×3 Matrix{Float64}:
 0.0182632   0.0550691  0.0995698
 0.00443534  0.0847981  0.0302305
 0.0138278   0.0285697  0.130521posterior mean  
3×3 Matrix{Float64}:
 -1.13067    1.064      2.19547
  0.274592  -1.62134    0.672081
  0.856081   0.557335  -2.86755

We can now sample from the BayesianGenerator object which gives a random matrix drawn from the posterior distribution

rand(Q_bayes)

3×3 Matrix{Float64}:
 -1.28443    0.987111   1.46499
  0.308094  -1.50439    0.778556
  0.976332   0.517283  -2.24355

If we call the rand function again we get a different realization

rand(Q_bayes)

3×3 Matrix{Float64}:
 -1.09693    1.25391   1.9438
  0.239546  -1.93427   0.746049
  0.857388   0.68037  -2.68985

We can call the rand function with an additional integer argument to get a list of realizations

number_of_samples = 100
Q_bayes_list = rand(Q_bayes, number_of_samples)

100-element Vector{Matrix{Float64}}:
 [-1.0863938382287588 1.090269012912725 1.8294736565302356; 0.19034429494915964 -1.4501370816881989 0.7514547478446989; 0.8960495432795991 0.3598680687754736 -2.5809284043749345]
 [-1.0481184776034607 0.9615557733566846 1.5616337274921526; 0.18828311181296997 -1.4565290107987152 0.8950535312730717; 0.8598353657904908 0.49497323744203053 -2.4566872587652244]
 [-1.2698598105757197 1.2208312736856268 1.937249429099191; 0.20054011420791645 -1.8903568361757102 0.8759461106502346; 1.0693196963678033 0.6695255624900834 -2.8131955397494255]
 [-1.1888853760777647 1.010365312213905 2.2146059784474006; 0.39523536938456827 -1.2915796888818538 0.42681276687190384; 0.7936500066931963 0.281214376667949 -2.6414187453193043]
 [-1.0197393352742763 1.256930281016118 2.169516964930346; 0.1900439358417498 -1.621802548515956 0.6702831893882244; 0.8296953994325263 0.3648722674998379 -2.8398001543185707]
 [-1.3566801253059724 0.9836714876731875 1.9971110907260377; 0.4544932732165382 -1.9954717125432224 0.6411853964437394; 0.9021868520894342 1.0118002248700346 -2.638296487169777]
 [-0.9296887754457219 1.2319346163589049 2.009565064606964; 0.23271405802556192 -1.673018972597842 0.8424160525918357; 0.69697471742016 0.4410843562389371 -2.8519811171988]
 [-0.968913671317772 1.0896148711124527 2.144405558458069; 0.2367978009922009 -1.4241143540854218 0.6140979792533799; 0.732115870325571 0.3344994829729689 -2.758503537711449]
 [-1.0340242045139374 0.8913971102385871 2.430519234716; 0.13644840607754763 -1.4109267726784902 0.6508517999410508; 0.8975757984363898 0.5195296624399028 -3.081371034657051]
 [-1.4461887692959006 1.4819398253428995 2.0072604442776716; 0.3588824217190628 -2.4305308494575555 0.7069475287375584; 1.0873063475768376 0.9485910241146561 -2.7142079730152298]
 ⋮
 [-0.9905349304660338 0.7903643365852441 3.042914552242634; 0.24956093628379364 -1.5915147218019678 0.5422118466868081; 0.7409739941822401 0.8011503852167237 -3.5851263989294417]
 [-1.4380052624532704 1.236037795853921 1.8225545016889548; 0.34058681266369917 -1.8479577949377066 1.0824841704295973; 1.0974184497895714 0.6119199990837852 -2.905038672118552]
 [-1.3487498703705567 0.7411020111619587 2.747152004778648; 0.2946371639280493 -1.13471466422307 0.8373995151017881; 1.0541127064425075 0.3936126530611113 -3.5845515198804367]
 [-1.0429776786882903 0.9725809131522108 2.687974888295159; 0.18617471730961602 -1.6030116413987188 0.5510271710843626; 0.8568029613786743 0.6304307282465076 -3.2390020593795215]
 [-1.0424250152841803 1.6588470881442705 1.6906362692966033; 0.29053663813110514 -2.1984700800332986 0.4705070271405873; 0.7518883771530751 0.539622991889028 -2.1611432964371904]
 [-1.146297472605456 1.2036598098568059 2.504736422553042; 0.1572237263062368 -1.6321023546951519 0.8529428666926973; 0.989073746299219 0.42844254483834615 -3.35767928924574]
 [-1.055211772995256 0.8965720052753188 2.3607594828334686; 0.3270545104511646 -1.3402075920048198 0.40747508743530936; 0.7281572625440913 0.443635586729501 -2.768234570268778]
 [-1.001712564109727 1.1483245085694571 2.91762948822196; 0.23958667143228454 -1.5164472509992089 0.8978623632055582; 0.7621258926774426 0.36812274242975157 -3.8154918514275185]
 [-1.1909510261909961 0.9165353994206388 2.5309066496901798; 0.24645963883016914 -1.316052336119224 0.7090360612213963; 0.9444913873608269 0.3995169366985852 -3.2399427109115755]

from whence we can compute the sample mean and compare to the analytic mean

mean(Q_bayes_list) - mean(Q_bayes)

3×3 Matrix{Float64}:
  0.00725497   0.0339316   0.0183639
 -0.00319942  -0.0197629   0.0236872
 -0.00405555  -0.0141687  -0.0420511

and similarly for the sample variance

var(Q_bayes_list) - var(Q_bayes)

3×3 Matrix{Float64}:
 0.00230174    0.00298935  0.0192147
 0.000521285   0.00155187  0.00152409
 0.000228445  -0.00171988  0.00283814